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ABSTRACT 

We describe planetesimal accretion calculations in the Kuiper Belt. Our 
evolution code simulates planetesimal growth in a single annulus and includes 
velocity evolution but not fragmentation. Test results match analytic solutions 
and duplicate previous simulations at 1 AU. 

In the Kuiper Belt, simulations without velocity evolution produce a single 
runaway body with a radius ^ 1000 km on a time scale oc Mq^Cq, 
where Mq is the initial mass in the annulus, eo is the initial eccentricity of the 
planetesimals, and x ^ 1-2. Runaway growth occurs in 100 Myr for Mq ^ 10 
M0 and eg ~ 10^^ in a 6 AU annulus centered at 35 AU. This mass is close 
to the amount of dusty material expected in a minimum mass solar nebula 
extrapolated into the Kuiper Belt. 

Simulations with velocity evolution produce runaway growth on a wide range 
of time scales. Dynamical friction and viscous stirring increase particle velocities 
in models with large (8 km radius) initial bodies. This velocity increase delays 
runaway growth by a factor of two compared to models without velocity 
evolution. In contrast, collisional damping dominates over dynamical friction 
and viscous stirring in models with small (80-800 m radius) initial bodies. 
Collisional damping decreases the time scale to runaway growth by factors of 
4-10 relative to constant velocity calculations. Simulations with minimum mass 
solar nebulae, Mq ~ 10 Mq, and small eccentricities, e ^ 10~^, reach runaway 
growth on time scales of 20-40 Myr with 80 m initial bodies, 50-100 Myr with 
800 m bodies, and 75-250 Myr for 8 km initial bodies. These growth times vary 
linearly with the mass of the annulus, oc Mo"\ but are less sensitive to the 
initial eccentricity than constant velocity models. 

In both sets of models, the time scales to produce 1000+ km objects are 
comparable to estimated formation time scales for Neptune. Thus, Pluto-sized 
objects can form in the outer solar system in parallel with the condensation of 
the outermost large planets. 
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1. INTRODUCTION 



Current models for planet formation involve aggregation of solid planetesimals and 
gas accretion in a circumstellar disk (see, for example, [Hayashi et al. 1985| ; Poss 1993| ; 



and references therein). Large dust grains within the disk first settle to the midplane. 
These grains may then coagulate into successively larger grains (e.g., [Weidenschilling 1980| ; 
Weidenschilling fc Cuzzi 1985|) or continue to settle in a very thin layer that eventually 



becomes gravitationally unstable (e.g., [Goldreich fc Ward 1973| ). Both paths produce 



km-sized planetesimals that collide and merge to produce large bodies such as planets. 
Despite the complex and sometimes unknown physics, many simulations produce objects 
resembling known planets on time scales roughly comparable to the expected lifetime of the 



protosolar nebula (see, for example, |Safronov 1969|; |Greenberg et al. 1978|, 1984; 


IN akagawa 


et al. 19831; |Wetherill fc Stewart 1989 


; 1993; 


Spaute et al. 1991 


; Kolvoord fc Greenberp 


1992; 


Weidenschilling fc Davis 1992; Pollack et al. 1996). 



Recent observations of slow-moving objects in the outer solar system offer a new 
challenge to planet formation models. The trans-Neptunian region now contains several 
dozen Kuiper Belt objects (KBOs) with estimated radii of 100-300 km (|Jewitt et al. 1996| ). 
The orbits of known KBOs suggest a division into at least three dynamical components with 
an inner radius of 30 AU and an unknown outer radius: (i) the classical KBOs, objects with 
roughly circular orbits, (ii) the resonant KBOs, objects in orbital resonance with Neptune 
( [Jewitt et al. 1996| ), and (iii) the scattered KBOs, objects with large, eccentric orbits ( |Luu 



et al. 19971) . Although the known population is still small, Jewitt et al. (1996) estimate 
that the region between 30 and 50 AU contains ~ 70,000 objects larger than 100 km. 
The total mass in the classical Kuiper Belt is thus at least 0.1 Mq. This mass probably 
represents a small fraction of the initial mass, because dynamical interactions with Neptune 
reduce the number of KBOs on short time scales compared to the age of the solar system 
( ILevison fc Duncan 199'^ ; [Malhotra 19961) . 

Despite these new observations, the origin of KBOs remains uncertain. Edgeworth 
(1949) and Kuiper (1951) first suggested that the Kuiper Belt was a natural extension of 
the original solar nebula. Holman fc Wisdom (1993) later showed that small KBOs, once 
formed, can survive at 30-50 AU for times approaching the age of the solar system. More 
recent dynamical studies confirm this conclusion and explain the observed distribution of 
KBOs in a general way ( [Levison fc Duncan 1993| ). The formation process and time scale 
for KBOs, however, is still controversial. Planetesimal simulations for plausible protosolar 
nebulae at 25-30 AU show that Neptune can grow to its present size in 10-100 Myr 
( Fernandez fc Ip 1981|, 1984; |Ip 1989| ; [Pollack et al. 1996|) . These results suggest that Pluto 
might form on a similar time scale at ~ 40 AU, because growth times are not a steep 
function of semi-major axis (see [Aarseth et al. 1993 ; Pollack et al. 1996 ; and references 
therein). Nevertheless, Stern (1995, 1996) and Stern fc Colwell (1997) conclude that KBO 
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formation requires 100-1000 Myr for the conditions expected in the outer solar system. 
These hmits far exceed the time scale required to produce Neptune, whose accretion time 
is constrained by the 10-100 Myr lifetime of the protosolar nebula (see [Pollack et al. 1996 
and references therein). 

In this paper, we attempt to resolve the uncertainties surrounding KBO production 
with a new planetesimal simulation at 35 AU. We have developed an evolution code to 
follow the growth and velocity evolution of planetesimals with a wide range of initial masses. 
The code matches analytic models and duplicates Wetherill & Stewart's (1993; hereafter 
WS93) simulation of planetesimal evolution at 1 AU. Our numerical results demonstrate 
that small bodies with initial radii of 80 m to 8 km can produce 1000+ km objects on time 
scales of 10-100 Myr. We confirm these calculations with a simple analytic estimate of the 
growth time as a function of semi-major axis. This analysis supports previous estimates for 
a short growth phase for Neptune, 10-100 Myr, and indicates that Pluto-Charon can form 
just outside the current orbit of Neptune on a similar time scale. 

We outline the accretion model in Sec. 2, describe our calculations in Sec. 3, and 
conclude with a discussion and summary in Sec. 4. The Appendix contains a complete 
description of the algorithms and detailed comparisons with analytic models. 



2. The Accretion Model 

For our simulations of accretion in the Kuiper Belt, we adopt Safronov's (1969) 
particle-in-a-box method, where planetesimals are treated as a statistical ensemble of 
masses with a distribution of horizontal and vertical velocities about a Keplerian orbit. Our 
simulations begin with a differential mass distribution, n(mj), in a single accumulation zone 
centered at a heliocentric distance, a, with an inner radius at a — Aa/2 and an outer radius 
at a + Aa/2. We approximate the continuous distribution of particle masses with discrete 
batches having particle populations nj(t) and total masses Mj(t) (WS93). The average mass 
of a batch, mi{t) = Mi{t)/ni{t), changes with time as collisions add and remove bodies 
from the batch. This procedure naturally conserves mass and allows a coarser grid than 
simulations with fixed mass bins ( [Wetherill 1990| , and references therein; [WS93[) . 



To evolve the mass and velocity distributions in time, we solve the coagulation and 
energy conservation equations for an ensemble of objects with masses ranging from ~ 10^^ g 
to ~ 10^^ g. The Appendix describes our model in detail and compares our numerical 
results with analytic solutions for standard test cases. We adopt analytic cross-sections to 
derive collision rates and compute velocity changes from gas drag and collective interactions 
such as dynamical friction and viscous stirring. Our initial approach to this problem 
ignores fragmentation, which we will consider in a later paper. In this study, we focus on 
developing a good understanding of planetesimal growth as a function of initial conditions 
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in the Kuiper Belt. 

To test our numerical procedures in detail, we attempt to duplicate WS93's simulations 
of planetary embryo formation at 1 AU. WS93 (see also [Wetherill fc Stewart 1989| ; [Barge 
fc Pellat 1990| , 1991, 1993; ^paute et al. 1991| ; [Aarseth et al. 1993|) demonstrate that an 
ensemble of 8 km objects can produce 10^^ g (Moon-sized) objects on a 10^ yr time scale. 
WS93's model begins with 8.33 x 10^ planetesimals having radii of 8 km and a velocity 
dispersion of 4.7 m s~^ relative to a Keplerian orbit (Table 1; see also Table 1 of WS93 ). 
Tables 2-3 summarize our results using the WS93 initial conditions with a mass spacing 
factor oi 6 = mi^i/mi = 1.25 and 1.4 between successive mass batches and two different 
analytic cross-sections. Figure 1 shows our reproduction of the WS93 results without 
fragmentation for 6 = 1.25. This simulation produces 14 3-9 x 10^^ g objects in 1.5 x 10^ 
yr, which agrees with the results in WS93 (see their Figure 12). Our simulation confirms 
the broad "plateau" in the cumulative number, Nc, at log = 24-26 and the rough 
power law dependence, Nc oc m^^, at log = 21-23. The broad plateau extends across 
a smaller mass range and becomes more rounded as 6 increases (Figure 2). The maximum 
planetesimal mass, rrimax, at the conclusion of the calculation depends on both 6 and the 
cross-sections. We find marginally larger rrimax for the Spaute et al. (1991) cross-sections. 
In general, rrimax increases as S increases. 

The evolution of particle velocities in our simulations agrees with the WS93 results 
(Figure lb). All of the velocities increase monotonically with time due to viscous stirring. 
The velocities of the larger bodies increase very slowly, because dynamical friction transfers 
their kinetic energy to the smaller bodies. The simulation maintains a nearly constant ratio 
of vertical to horizontal velocity, Vi/hi ~ 0.53, for all but the most massive bodies, which 
have Vi/hi < 0.5. The equilibrium ratio of Vi/hi yields < i > / <e>~0.6, in agreement 
with Barge & Pellat (1990, 1991; see also Pornung et al. 1985| ). At the conclusion of the 
simulation, our velocities for small bodies, hi ^ 500 m s~^ at m, ~ 10^^ g, are roughly 50% 
larger than those of WS93. Our velocities for large bodies, /ij ^ 10 m s~^ at rrii ~ 10^^ g, 
are roughly 50% smaller than those of WS93. We also fail to reproduce WS93's abrupt 
drop in hi at log rui = 24. However, these differences - which are independent of ^ - have 
a negligible effect on the final mass distribution and probably result from slightly different 
algorithms for low velocity collisions. 

Gas drag is included in our simulations but has a negligible impact on the evolution. 
All of the 1 AU models lose ~ 0.01% of their initial mass over 1.5 x 10^ yr. Velocity 
changes due to gas drag are essentially zero, because the particle masses are so large. 

To understand the sensitivity of these results to initial conditions, we consider 
the growth time of planetesimals from the coagulation equation, equation A4. 
For most cases of interest, the growth time for bodies with rUi is approximately 
r ^ nQ/{dn/dt) oc HaAa/njVFg{ri + rj)^, where rij is the number of lower mass bodies. 
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V is the relative velocity, Fg is the gravitational focusing factor, rj and vj are the radii of 
particles i and j, and H is the vertical scale height. Collisions between low mass objects 
are in the high velocity regime, where the gravitational focusing factor is F„ ^ 1 and 



T oc a^/^ Aa 



(n + r,- 



This growth time is independent of the initial e and i. 



Gravitational focusing becomes effective in low-velocity collisions of massive objects; the 
growth time then depends on the initial velocity and is r oc a^/^ Aa nj^ (r^ + rj)~^ V"^. 
The extreme sensitivity of the growth time to velocity is the reason why low-velocity 
planetesimals experience runaway growth in our 1 AU simulations ( |Wetherill fc Stewart 



|l989t [Ida fc Makino 1992a| , 1992b; [Kokubo fc Ida 1996| ; see also [WS9^ and references 
therein). We adopt 1000 km as a useful reference radius and write the time for 8 km objects 
to produce 1000 km objects at 1 AU as (see also [Barge fc Pellat 1990| ): 



^0 



a 

1 AuJ 



Aa 



Po 



1/3 



0.17 AU/ 13 g cm- 



450 cm s 



6 X 10^ ^g 

, Mo 



Using our simulations with 6 = 1.4, we derive the proportionality constant for this standard 
case with velocity evolution, r^^y = 10700 yr, and for a model with no velocity evolution, 
To,nv = 3750 yr. Additional simulations confirm the mass, velocity, and volume dependence 
of this relation for factor of two variations in a, 6a, po, Vq, and Mq about the values in Eqn 
1 (see also jAarseth et al. 19"93[; Pollack et al. 1996|). 



3. Kuiper Belt Calculations 



3.1. Starting Conditions 



To choose appropriate constraints on planetesimal simulations in the outer solar 
system, we rely on observations of other stellar systems and models of the protosolar nebula. 
First, current data indicate lifetimes of ~ 5-10 Myr for typical gaseous disks surrounding 
nearby pre-main sequence stars and for the solar nebula ([Sargent fc Beckwith 199"3[ ; [StromI 



\et al. 1993| ; [Russell et al. 19"9^ ). We adopt this estimate as a rough lower limit to the 



formation time scale of KBOs and assume that interactions between gas and planetesimals 
disappear on a similar time scale, Tg 10 Myr (see the Appendix, equation A29). Neptune 
formation places an upper limit on the KBO growth time, because Neptune excites KBOs 
through gravitational perturbations. Recent calculations suggest Neptune can form in 
5-100 Myr (Ip 1989; Lissauer & Stewart 1993; Lissauer et al. 1995; Lissauer et al. 1996). 
Once formed, Neptune inhibits KBO formation at 30-40 AU by increasing particle random 
velocities on time scales of 20-100 Myr ( [Holman fc Wisdom 1993| ; Puncan et al. 19"90[ ). We 



thus adopt 100 Myr as a rough upper limit to the KBO formation time scale at 30-40 AU. 
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We assume a wide range of starting conditions for KBO simulations. Our model 
annulus is centered outside the orbit of Neptune at 35 AU and has a width of 6 AU. This 
annulus can accommodate at least 10-100 isolated bodiesQ with nii ^ lO^'* g for e < 0.01. 
The simulations begin with Nq bodies of radius vq, with tq = 80 m, 800 m, and 8 km. 
These bodies have small initial eccentricities, e ~ 10"'^ to 10^^ ( |Malhotra 1995| ), and an 
equilibrium ratio of inclination to eccentricity, f3 =< i > / < e > = 0.6 ( Parge fc Pellat 



|1990| , 1991, 1993). The mass density of each body is fixed at 1.5 g cm~^. To set Nq, we 
extend the minimum mass solar nebula to the Kuiper Belt and integrate the surface density 
distribution for solid particles, S = So(a/ao)~^^^) across the 6 AU annulus. The dust mass 
is then Mmin ~ 0.25 Sq at 32-38 AU (ao = 1 AU). Most minimum mass solar nebula 
models have Sq = 30-60 g cm~^, which sets Mmm ~ 7-15 ( |Weidenschilling 1977 ; 
Hayashi 1981| ; [Bailey 1994| ). We thus consider models with initial masses of Mq = 1-100 



M0 to allow for additional uncertainty in Eq. Table 1 compares input parameters for all 
Kuiper Belt models with initial conditions at 1 AU (see also |WS93| ). Tables 4-5 summarize 
other initial conditions and results for the Kuiper Belt simulations summarized below. 

Our success criteria are based on direct observations of KBOs. The present day Kuiper 
Belt contains at least 70,000 objects with diameters exceeding 100 km at 30-50 AU ( |Jewitt| 
|fc Luu 1995| , [Jewitt et al. 1996] ). This population is some fraction of the initial Kuiper Belt 
population, because Neptune has eroded the Kuiper Belt over time (Polman fc Wisdom 
|1993|; puncan et al. 1990|) . Thus, a successful KBO simulation must achieve ^ 50 km 
in ^ 100 Myr, where is the radius where the cumulative number of objects exceeds 
Nc > 10^. Pluto formation is our second success criterion: plausible models must produce 1 
or more objects with maximum radius Tmax > 1000 km. In models with velocity evolution, 
we end simulations at 100-200 Myr or when r^aa; exceeds ~ 1000 km. To evaluate the 
dependence of runaway growth on the initial conditions, we extend simulations without 
velocity evolution to 5000 Myr or to when r-max exceeds ~ 2000 km. 



3.2. Models Without Velocity Evolution 

To isolate important processes in trans-Neptunian planetesimal evolution, we begin 
with constant velocity solutions to the coagulation equation. We ignore fragmentation and 
fix the velocities for all masses at hi = 4.0 emit km s~^ and Vi = 3.6 sin iinn km s~^. The 
initial eccentricity and inclination are set at iinn = 0.60 einn- The total mass and kinetic 
energy remain constant throughout the calculation. We also adopt a coarse mass spacing 



^"Isolated bodies" are planetesimals that cannot collide with one another, as defined in 
the Appendix following Eqn (A5) 
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factor, S = 1.4. This choice hmits our abihty to follow runaway growth with high accuracy 
during the late stages of the simulation but allows us to investigate a wide range of initial 
masses and velocities with a modest investment of computer time. Finally, we adopt simple 
formulae for gravitational focusing to speed our calculations (see equation A18; ppaute 
al. 19911) . Table 4 summarizes the initial conditions and results for models with Mq = 

init 



1-100 Mm and e,„H = 10"^ and 10 



Figure 3 shows how Nq evolves with time for a model with an initial planetesimal 
radius ro = 8 km, total mass Mq = 10 Mq, and Cinit = 10^^. This model begins with 
1.87 X 10^° initial bodies and produces ~ 42,500 objects with twice the initial mass after 
r = 100 yr. Roughly half of the original population experiences at least one collision by 
r ^ 16 Myr, when the 17 largest bodies have rrij ^ 10^° g. Slow growth continues until 
r ^ 59 Myr when the 3 largest objects have sizes comparable to large KBOs, rj ^ 100 km 
and rrii ^ 10^^ g. The growth rate of the large masses then increases considerably due to 
gravitational focusing. Runaway growth ensues. The cumulative mass distribution then 
follows a power law, Nc oc at low masses and develops a high mass shoulder that 

extends to larger and larger masses as the simulation proceeds. This shoulder resembles 
the runaway plateau observed in 1 AU models but does not evolve into a true plateau 
with Nc ~ constant as in Figs. 1-2. The largest planetesimals reach Vmax ~ 200 km at 
T ^ 66 Myr; Vmax exceeds 1000 km only 9 Myr later. A single runaway body with Vmax ~ 
4000 km begins to sweep up lower mass planetesimals at r ^ 80 Myr; by r ~ 85 Myr, it 
contains essentially all of the mass in the annulus. 

Simulations with tq = 8 km produce runaway growth independent of the initial mass in 
the annulus. Figure 4(a) indicates that each model experiences a long, linear growth phase 
until r^nax ~ 100-200 km. The largest objects then begin a short, rapid growth phase that 
produces several isolated, runaway bodies with Vmax ~ 1000 km. These runaway bodies 
accumulate all of the lower mass bodies and may merge to form a single runaway body if 
the isolation criterion permits. The time to produce runaway bodies with rj = 1000 km 
scales with the mass in the annulus, Tr ~ 753 Myr (Mq/I Mq)"^. For comparison, our 
scaling relation in Eqn (1) predicts Tr ~ 775 Myr (Mq/I M0)~^ for po = 1-5 g cm~^ in a 6 
AU annulus centered at 35 AU. 

Runaway growth also occurs independent of the initial radius, tq (Figures 4(a)-4(c)). 
Due to smaller initial cross-sections, models with ro = 80-800 m take longer to reach the 
rapid growth phase. These models make the transition from rapid growth to runaway 
growth more quickly, because it is easier for 100+ km objects to sweep up small objects with 
r < 1 km. In all cases, a single runaway body with r > 1000 km eventually accumulates all 
of the mass in each simulation, although the time scale is quite long, ^ 2700 Myr, for 
simulations with Mq = 1 M^ and ro = 80 m. Again, the runaway growth time scales with 
mass: Tr ^ 2340 Myr (Mo/1 M^)-^ for ro = 800 m and Tr ^ 2700 Myr (Mq/I M^)'^ for 
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tq = 80 111. The small increase in with initial radius for tq ^ 800 m suggests that models 
with ro < 80 m will reach runaway growth on time scales of ~ 3000 Myr, which is ~ 40 
times slower than models with ro = 8 km. 

Our results also confirm the velocity dependence derived in Eqn (1). Low eccentricity 
simulations with 50% smaller initial velocities reach runaway growth in 25% of the time for 
our standard model; simulations with twice the initial velocity require four times as long 
to achieve runaway growth. This simple relation begins to break down as the eccentricity 
increases to e ~ 10~^, as outlined below. The runaway time also scales with the width of 
the annulus, Aa, and the semi-major axis, a, as indicated in equation (1). 

High eccentricity models also achieve runaway growth but do not follow precisely the 
velocity scaling in Eqn (1). Figure 5 shows the radius evolution for models with various 
ro and ttiq for e = 10~^ (see also Table 4). The growth time for 1000+ km objects is 
Tr ~ 20-25 Gyr (Mo/1 M^)"^ nearly independent of the initial radius and velocity. This 
relation contrasts with the low eccentricity results, where the growth time is very sensitive 
to the initial conditions. In all of our simulations, planetesimal growth is orderly until 
gravitational focusing becomes important and runaway growth occurs. However, the radius 
where gravitational focusing becomes important increases from rj ~ 10 km at e = 10~^ 
to rj ~ 100 km at e = 10~^. For models with small initial bodies, ro ^ 800 m, the time 
scale to reach runaway growth is directly proportional to e. For models with large initial 
eccentricity, the long orderly phase also "erases" memory of the initial radius. Thus, t,. 
is nearly independent of ro for large e. The relatively short orderly growth phase of low 
e models does not erase memory of ro; Tr decreases with increasing ro for Cinu ^ 0.05. 
For models with ro ~ 8 km, gravitational focusing accelerates growth immediately at low 
e. These simulations do not have an orderly growth phase; instead, they follow the 1 AU 
simulations and satisfy the scaling relation in Eqn (1). 

Before we consider Kuiper Belt simulations with velocity evolution, our basic result 
that constant velocity models achieve runaway growth deserves some comment. First, 
previous simulations at 1 AU show that runaway growth requires dynamical friction to 
decrease the velocities of the largest bodies to the regime where gravitational focusing 
becomes important ( |Wetherill fc Stewart 1989| , 1993; |Barge fc Pellat 1990| , 1991; |lda 1990 



Ida fc Makino 1992a| , 1992b; [Ohtsuki 199^ ; |Kokubo fc Ida 1996| ). In Kuiper Beh models 
with einit = 10~^, gravitational focusing factors become very large at planetesimal masses 
of 10^^ to 10^^ g. Further growth of these bodies only enhances gravitational focusing, 
because the escape velocity increases while the impact velocities remain low. More massive 
objects thus "run away" from their lower mass counterparts. This response occurs in any 
constant velocity simulation as long as bodies can reach masses where the escape velocity is 
large compared to the relative impact velocity, Ve^ij/Vij 3> 1. Models with Cinu ^0.1 never 
reach this limit for plausible Mq and thus do not experience runaway growth. Models with 
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^init ~ 0.05 always produce runaway bodies, albeit at much later stages than models with 
Cinit ~ 10"^ (Figure 5). 

Our final comment on runaway growth concerns the shape of the cumulative number 
distribution near the end of the simulation. During runaway growth, models at 1 AU 
develop a plateau in the cumulative number distribution that extends from rrii = 10^^ g 
to rui = 10^^ to 10^^ g (see Figures 1-2). This plateau separates runaway bodies from 
the lower mass objects which are in the orderly growth regime and have a power law size 
distribution, Nc oc r^^ for log = 21-24 ( |WS93| ; see Figures 1-2). The Kuiper Belt 
simulations also produce a power law size distribution, Nc oc r^'^''^ for rrii ^ 10^^ g, but 
they develop a high mass "shoulder" instead of a marked plateau at runaway growth (see 
Figure 3). To test if this feature is a function of the mass resolution as in 1 AU models, we 
simulate evolution at 35 AU with 5=1.1 and 1.25 for Mq = 10 and tq = 8 km. As the 
mass resolution in the simulation increases from 6 = 1.4 to 6 = 1.1, the high mass shoulder 
follows a very shallow power law, Nc oc r"^''' (Figure 6). This power law becomes better 
defined as the mass resolution increases further, but it never develops into the "runaway 
plateau" produced in the 1 AU models (Figures 1-2). This result suggests that the broad 
plateau in 1 AU models is the result of velocity evolution, which reduces the velocity of the 
most massive objects and enhances gravitational focusing ( |WS93|) . We will now test this 
hypothesis by considering Kuiper Belt models with velocity evolution. 



3.3. Models with Velocity Evolution 

To understand the importance of velocity evolution in the Kuiper Belt, we add several 
physical processes to the calculation: (i) gas drag, (ii) dynamical friction and viscous 
stirring from long-range (elastic) collisions, and (iii) dynamical friction and viscous stirring 
from short-range (inelastic) collisions. As in our constant velocity models, we begin with A^o 
bodies at radii, ro = 80 m, 800 m, and 8 km. We adopt 6 = 1.4 for the mass spacing factor 
and use our simple expression for gravitational focusing, Eqn A18. The initial velocities 
are hi = 4.0 (ej„jt/10"^) m and Vi = 2.1 (ej„ji/10~^) m s~^. The eccentricity and 
inclination evolve separately due to collisions and collective interactions (see section A. 3 of 
the Appendix). Table 5 summarizes the initial conditions and results for models with Mq 
= 1-30 M0 and Cinu = lO'^ to lO'^. 

Before describing the results of our simulations, it is useful to compare various time 
scales for velocity evolution at 35 AU. First, gas drag is negligible in models that ignore 
fragmentation. A typical simulation at 35 AU loses ~ 10^^% of its total mass due to gas 
drag in 100 Myr. Velocity changes due to gas drag are also insignificant; the time scale for 
gas drag to modify the velocity exceeds 10 Gyr for all masses in our simulation. 

Velocity changes due to elastic and inelastic collisions, however, are significant. Figure 
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7(a) compares time scales, t„_/j = hi/{dhi/dt), for horizontal velocity evolution as a function 
of particle mass at 35 AU. The two curves show r^^h for interactions between particles 
of the same mass with a power-law size distribution, Nc oc r"^''^, and constant velocity. 
These time scales are not integrated over the size distribution and are relevant only when a 
simulation has a small range of masses. Viscous stirring - which tends to increase particle 
velocities - is ineffective for rrii ^ 10^^ kg, r^^h ~ 10^^ yr, but very effective, r^^h ^ 10^ yr, 
at rrii ^ 10^'^ kg. CoUisional damping is also more effective at large masses, but the time 
scale is much less mass-sensitive than viscous stirring. CoUisional damping balances viscous 
stirring for an initial particle mass, mo ~ 10^^ kg, which corresponds to Tq ~ 5 km. 

To illustrate these points in more detail. Figures 7(b) and 7(c) plot the integrated time 
scales, T^^h = J2f=i hi / {dhi / dt) , for the horizontal velocity at two stages of a model with 
velocity evolution. In Figure 7(b), the maximum mass has m^ax ~ 10^^ kg. CoUisional 
damping still dominates viscous stirring for the lowest masses, but the mass where the two 
processes balance has moved from rui ~ 10^^ kg to ~ 5 x 10^^ kg. Once rumax ~ 10^^ 
kg, viscous stirring dominates coUisional damping for all masses. Particle velocities thus 
increase once massive objects with ^ 100 km are produced. The time scale for viscous 
stirring is quite short, ^ 10^ yr, at the lowest masses considered in our models, so the 
velocity increases can be large during the 100 Myr time scale of a typical simulation. 

Figure 8 shows how Nc and hi evolve with time in a model with tq = 8 km, Mq — 
10 M0, and Cinit = 10^'^. The simulation begins with Nq = 1.87 x 10"'^'^ and produces 
762 objects having eight times the initial mass in 1 Myr. Roughly half of the initial bodies 
experience at least one collision by 24 Myr, when the largest object has r^^x ~ 37 km. The 
horizontal velocities then range from hi — 1 m s~^ at n = 37 km up to /ij = 7 m s~^ at 
ri — 8 km. Orderly growth produces = 100 km objects at 56 Myr, and this population 
reaches Nc ~ 100 at 61 Myr. This phase continues until ~ 100 Myr, when Charon-sized 
objects with ^ 500 km begin to grow rapidly. There arc 10 "Charons" at 110 Myr, 47 
"Charons" at 125 Myr, 107 "Charons" at 150 Myr, and 202 "Charons" at 180 Myr when we 
ended the simulation. At 180 Myr, 20 Pluto-sized objects with ~ 1000 km are isolated 
bodies about to run away from the rest of the mass distribution. 

In contrast to the constant velocity simulations with low Cinu, this model does not 
immediately enter a rapid growth phase once objects with ~ 100 km are first produced. 
Viscous stirring and dynamical friction slowly increase the velocities of low mass bodies 
throughout the simulation: the horizontal velocity increases from hinn = 4 m s~^ to hi ~ 
65 m s~^ at 180 Myr (Figure 8). This twcntyfold increase in the eccentricity reduces 
gravitational focusing by a factor of 400 and retards the growth of the most massive objects. 
Evolution thus proceeds at a pace between the constant velocity simulations with einit — 
IQ-^' and Cinu = lO'^. 

The evolution for simulations with smaller initial masses is different, because coUisional 
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damping then dominates the velocity evolution (Figure 7). Figure 9 shows the time 
evolution of Nc and hi for ro = 800 m, Mq = 10 Mq, and Cinit = 10~^. This simulation 
begins with Nq = 1.87 x 10^^ and produces five objects having eight times the initial 
mass in 1 Myr. It takes only 7.7 Myr for half of the initial objects to collide at least 
once. The maximum radius is then Vmax ~ 2.5 km. Velocity damping from inelastic 
collisions overcomes viscous stirring, so the particle velocities remain low and do not change 
significantly with mass. As evolution proceeds, dynamical friction efficiently damps the 
velocities of the largest bodies, but coUisional damping still maintains modest velocities at 
low masses. Bodies with = 8-10 km begin to form at 32-33 Myr, when only 11% of the 
original objects remain. The evolution soon overtakes the tq = 8 km model. Orderly growth 
produces 50 km objects at 45 Myr and 100 km objects at 48 Myr. Runaway growth begins 
shortly thereafter. Charon-sized objects form at 60 Myr and reach Pluto-size at 80-81 Myr. 

Simulations with ro = 80 m reach runaway growth on even faster time scales. Figure 
10 shows the time evolution of Nc and hi for Mq = 10 M^ and Cinu = 10~^. At 1 Myr, only 
43% of the initial bodies have yet to experience a coUision; 33 objects already have 
270 m. The maximum radius reaches Vmax — 800 m in 9 Myr and Vmax — 8 km in 17 Myr. 
The low mass bodies first lose ~ 50% of their initial velocity, hinit = 4 m s~^, and begin to 
increase in velocity at 17-18 Myr when viscous stirring from long-range collisions finally 
overcomes damping from inelastic collisions. At this time, the high mass bodies have low 
velocities due to dynamical friction, hi ~ 0.01 m s~^, and begin to grow rapidly. A runaway 
plateau in the Nc distribution develops at 24 Myr and extends to Charon-sized objects at 
25 Myr. At the conclusion of this simulation at 33 Myr, five Pluto-sized objects have n 
900-1000 km. The velocities of these high mass objects are then hi ^ 0.03-0.05 m s~^. 

With their long runaway growth times, models with ro — 8 km cannot meet both of our 
success criteria unless the initial mass is very large, Mq ^ 15-20 M^. Viscous stirring and 
dynamical friction increase the velocities of the small objects throughout these simulations, 
which reduces gravitational focusing and delays runaway growth compared to models with 
smaller tq. The long approach to runaway growth allows the production of many large 
KBO's; simulations with Mq ~ 6-20 Mq have rs ~ 50-100 km and thus reach our first 
success criterion. However, the time scale to produce 1000-1- km objects is 4-5 times longer 
than models with ro = 80 m, i.e., fa 130 Myr (Mq/IO Mq)-^ Most of these models thus 
fail to make Pluto on a reasonable time scale. 

Models starting with lower mass objects, ro = 80 m and 800 m, meet both success 
criteria. Although viscous stirring and dynamical friction stir up the velocities of the lowest 
mass objects, the time scale for the velocity to increase is large compared to models with ro 
= 8 km (see Figure 7). The ro = 800 m models reach runaway growth faster than models 
with ro = 8 km and produce Pluto-sized objects in ~ 83 Myr (Mo/10 Mq)"^ for Cinu = 
10"^. The combination of smaller particle velocities and a shorter time to runaway growth 
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results in fewer KBOs compared to models with — 8 km. Nevertheless, these models 
achieve ~ 50-90 km during the runaway growth phase. 

Models with tq = 80 m and Cinu = 10~^ have the shortest runaway growth times and 
produce the fewest numbers of KBOs. The time scale to produce Pluto-sized planets is 
Tr ~ 32 Myr (Mq/IO Mq)^-*^ for Cinu = 10~^, which easily allows Pluto formation in a 
minimum mass solar nebula as Neptune forms at a smaller semi-major axis. These models, 
however, struggle to build a population of 10^ KBOs during the runaway growth phase. 
With relatively low particle velocities at all masses (see Figure 10b), objects with 
10-20 km do not grow as rapidly as larger bodies. This evolution tends to concentrate 
material in the more massive KBOs and reduces the number of lower mass KBOs with ~ 
50-100 km. This smaller in the calculations leads to partially-successful models that 
yield several Pluto-sized objects and ^ 50 km (Table 5). 

At large initial eccentricity, planctesimal growth follows the evolution of low eccentricity 
models but on longer time scales (Figure 12). In simulations with tq = 80 m and 800 
m, there are enough small bodies for inelastic collisions to damp the particle velocities 
substantially. Dynamical friction further decreases the velocities of the most massive bodies 
and allows runaway growth to occur on reasonable time scales. In simulations with ro = 80 
m and Cinu — 10~^, runaway growth is delayed by a factor of ~ 2.5 compared to simulations 
with Cinit — 10~^. This delay increases to a factor of ~ 3 for tq = 800 m. At tq = 8 km, 
collisional damping initially reduces particle velocities but is overcome by viscous stirring 
when r„iax ~ 50 km. At this time, the velocities arc large, ~ 30 m s~^, and growth 
is slow. None of these models reach runaway growth on a 100 Myr time scale. Runaway 
growth, if it occurs, is delayed by a factor of ~ 6 in large emit models compared to that in 
low Cinit models. 

Our results for large e^nit thus favor low mass initial bodies. Simulations with tq 
= 80 m and Cinit — 10~^ produce more KBOs with ~ 50-100 km than their low Cinit 
counterparts. The longer orderly growth phase and somewhat larger particle velocities at 
the onset of runaway growth favor the growth of 50-100 km objects. The for these mass 
distributions is 10% to 20% larger than the for low Cinit {r^ ~ 50-90 km for einn = 10^^ 
compared to ~ 50-75 km for e^nit = 10~^). Models with larger tq are less successful. For 
To = 800 m, runaway growth begins well after 100 Myr unless Mq ^ 20-30 M0. Runaway 
growth always occurs on a very long time scale, ^ 150-200 Myr, for ro — 8 km and Mq ^ 
30 M0. 

Unlike the constant velocity models, simulations with velocity evolution begin to 
develop a "runaway plateau" in the cumulative mass distribution when rmax ~ 500 km (see 
Figures 8-10). In constant velocity simulations, we found two power law cumulative mass 
distributions, A'',;^ oc r"^''''^ at low masses and A",^ oc r"-*^'^^ at large masses. The total mass 
per mass batch is then roughly constant at low masses and increases slowly with mass at 
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large masses. In models with velocity evolution, dynamical friction reduces the velocities 
of the largest bodies to ^ 0.1 m s~^ and maintains velocities of ~ 1-10 m s~^ for 
smaller objects with r, ^ 10 km. As noted by WS93, the increase in the escape velocity 
with mass coupled with the decrease in Vi produce substantial increases in the collisional 
cross-sections (see also [Wetherill fc Stewart 198^ ; [Barge fc Pellat 199q , 1991, 1993; |lda &| 
Makino 1992a| , 1992b, phtsuki 199^ ; [Kokubo fc Ida 1996|) . In our models, the velocity 
distribution resembles a step function and produces a step-like increase in the gravitational 
focusing factors, from Fg ~ 10-100 at n ~ 10-100 km up to Fg ~ 10^ at ~ 300-1000 
km. Runaway growth then converts the Nq oc r~^'^^ mass distribution into A^^ ~ constant, 
because objects with rj ~ 100-200 km grow too slowly to fill in the power law as objects 
with Tj ^ 500 km run away. At low masses, the size distribution remains a power law, 
Nc (X r~^'^^, because runaway growth does not change the size distribution significantly. 

To conclude this section. Figure 13 summarizes results for accretional evolution in 
the Kuiper Belt with velocity evolution and no fragmentation. Successful simulations that 
produce ~ 10^ KBOs and a few Pluto-sized objects on time scales of 100 Myr or less 
have initial masses somewhat larger than that predicted for a minimum mass solar nebula 
extrapolated into the Kuiper Belt, Mq ^ 10 Mq, and bodies with initial radii of tq ~ 
80 m to 8 km. Simulations with small initial bodies, tq = 80 m, tend to produce fewer 
KBOs than simulations with larger initial bodies, ro = 800 m and 8 km. This feature of 
our calculations results in partially-successful models that produce Pluto-sized objects but 
too few KBOs in 100 Myr. In models with Mq ?a 6 M0and Tq = 80 m, runaway growth 
removes KBOs from the mass distribution more rapidly than they are produced from lower 
mass objects. This evolution does not occur in models with larger initial bodies, because 
collisional damping is less effective at "circularizing" the orbits during the orderly growth 
phase. The higher particle velocities in these models allows formation of more KBOs during 
the runaway growth of Charon-sized objects. 

Collisional evolution often fails to produce 100+ km objects on any useful time scale. 
Simulations with Mq ^ 1-6 produce neither Pluto-sized objects nor a substantial 
number of 100+ km KBOs in 100 Myr. Large initial eccentricities exacerbate this problem 
for models with tq ^ 800 m, because collisional damping cannot reduce the particle velocities 
before 100+ km objects form. These simulations can produce KBOs and Plutos on longer 
time scales, 100-1000 Myr, in systems where a Neptune-sized object does not constrain the 
formation time. Extrapolating our results to smaller initial masses, simulations with Mq ^ 
0.1 M0 fail to produce KBOs during the age of the solar system, ~ 5 Gyr. 
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3.4. Limitations of the Models 

Statistical simulations of planetesimal growth have well-documented approximations 
and uncertainties. The model assumes a homogeneous spatial distribution of planetesimals 
whose velocities are small compared to the orbital velocity. These assumptions are good 
during the early stages of planetesimal evolution. As planetesimals grow, dynamical friction 
can reduce the velocities of high mass objects below limits where the statistical approach 
is valid (Parge fc Pellat 19901) . Once this limit is reached, runaway growth produces a few 



large bodies that are not distributed homogeneously in space (|WS93| ; [Kokubo fc Ida 19961) 



These large bodies can then pump up the velocities of the smallest bodies on short time 
scales through viscous stirring (Figure 7). We end the simulations with velocity evolution 
during the runaway growth stage when the basic assumption of a homogeneous distribution 
of planetesimals begins to break down. The velocities of low mass bodies remain small 
compared to the Keplerian velocity, but the most massive objects often have velocities 
below the low velocity limit of the kinetic approximation. We will discuss this problem 
below. 

The remaining limitations of the statistical approach involve our implementation of 
standard algorithms. We adopt a single accumulation zone and thus cannot follow the 
evolution in semi- major axis of a planetesimal swarm (see [Spaute et al. 1991]) . We use a 



coarser grid than some simulations, but this choice has little impact on the results at 35 
AU. At 1 AU, the lag of simulations with 5 > 1.1 relative to a simulation with 6 = 1.1 
increases with increasing 6; we find a 12% lag for 6 = lA but only a 2-3% lag for 6 = 1.25. 
At 35 AU, the lag in runaway growth relative to a 5 = 1.1 model increases from 4%-5% 
for S = 1.4 to 10%-15% at 5 = 2. Our S = 1.4 simulations thus overestimate the runaway 
growth time only by 4-5% (see also Wetherill 1990 ; Kolvoord fc Greenberg 1992|) . This 



error is small compared to other uncertainties in the calculation. 

Our choice of the initial mass distribution has a modest impact on our results. We 
calculated the evolution of several size distributions with equal mass per mass batch for 
^min ~ ri ^ Tmax- Simulations with rmm ~ 100-1000 m and Vmax ~ a few km are nearly 
indistinguishable from simulations with a single starting radius, tq ^ r^ax- In these models, 
collisional damping effectively reduces all particle velocities as described above and allows 
runaway growth to occur. Simulations with large Vmax ~ 8 km are similar to those with a 
single starting mass, unless r^m is small. For Vmin ~ 800 m, collisional damping keeps the 
particle velocities small compared to models with a single starting mass. Runaway growth 
occurs in these models, but the time scale to runaway is sensitive to r^m- We plan to 
explore this sensitivity in more detail when we include fragmentation in the calculation. 

The most uncertain approximation in our calculations is the treatment of low velocity 
collisions. During the late stages of most simulations, the massive bodies have very low 
velocities and very large gravitational ranges. The velocities are often smaller than the Hill 
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velocity, Vh, which invahdates the basic assumptions for our velocity evolution calculations 
( IBarse fc Pellat 199q , 1991, 1993). Barge & Pellat (1990) and WS93 have developed 
different approximations to low velocity collisions based on Ida's (1990) n-body simulations 
(see also [Ohtsuki 1992|) . We note that these two approximations produce different velocities 
for high mass planetesimals, which we plan to examine in more detail when we include 
fragmentation in our calculations. The mass evolution of runaway bodies is not affected by 
our treatment of velocity evolution in the low velocity regime. 

The large gravitational ranges of the most massive bodies also invalidate the standard 
treatment of collisions. We use the WS93 prescription for isolating the largest bodies from 
collisions with one another and adopt the Greenberg et al. (1992) approach to low velocity 
collisions in the two-dimensional regime. Removing the isolation criteria allows the largest 
body to grow more rapidly than the isolated bodies and reduces the time scale to runaway 
growth by 5%-10%. Removing the Greenberg et al. (1992) two-dimensional cross-sections 
has no substantial effect on our results. 

Aside from fragmentation, we have included all important physical processes in 
planetesimal evolution. Our neglect of fragmentation, however, is a serious limitation. In 
previous simulations, fragmentation of relatively strong bodies with ^ 1-10 km produces 
a significant amount of cratering debris that can be accumulated later by runaway bodies 
(e.g., [WS93| ; Parge &: Pellat 1993| ). This process usually becomes important only in the 
late stages of calculations at 1 AU: it slows growth during early stages but speeds up 
runaway growth later in the evolution (Parge fc Pellat 1993|; |WS93|). However, collisions 



between very weak bodies can disrupt and thereby prevent any growth of icy planetesimals 
at modest velocities. The importance of fragmentation at a ~ 35 AU thus depends on the 
unknown strength of KBOs. 

We can estimate the importance of fragmentation in Kuiper Belt simulations using 
Barge & Pellat's (1993) results for a reasonable fragmentation model. They adopt the 
Housen et al. (1991) energy prescription for planetesimal disruption and derive collisional 
outcomes for several test cases. These results are most appropriate for rocky asteroids, but 
it is straightforward to scale them to the weaker, icy bodies that might exist at 35 AU. We 
consider the two cases recently adopted by Stern & Colwell (1997): strong, rocky KBOs 
and weak, icy KBOs. 

Fragmentation does not significantly change our main conclusions if KBOs are strong 
objects. According to Figure 5 of Barge & Pellat (1993), fragmentation modifies the growth 
of 10 km bodies only when e ^ Ccrit ~ 0.025. The critical eccentricity for fragmentation 
decreases to ecrit ~ 10"^ for 1 km objects and Ccrit ~ 2 x 10"'^ for 0.1 km objects 
(see also Fig. 1 of [Stern 1996| ). Our low e simulations never reach these critical values. 
Fragmentation is important in large e simulations, but most of these models do not produce 
KBOs on a reasonable time scale. 
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The growth of icy KBOs is probably very sensitive to the time scale for velocity 
evolution. We expect fragmentation to dominate the early evolution of all simulations 
considered above, because only objects with ^ 20-30 km can survive collisions and 
produce larger bodies (see also Stern & Colwell 1997). As the evolution proceeds, however, 
inelastic collisions should damp the velocities of bodies with rj ^ 1 km, while dynamical 
friction damps the velocities of the most massive objects (see Figure 7). These damping 
time scales are short compared to the viscous stirring time scales, so the particle velocities 
decrease on relatively short time scales, ~ 10 Myr. This damping is probably sufficient to 
allow the growth of 1-10 km objects on time scales similar to those found in our models 
without fragmentation. Smaller bodies may not grow unless dynamical friction is very 
efficient. We will explore this possibility in our second paper. 



4. DISCUSSION AND SUMMARY 

We have developed a time-dependent planetesimal evolution program similar to the 
WS93 code used to simulate the formation of terrestrial embryos from small bodies. The 
program incorporates coagulation with realistic cross-sections and velocity evolution using 
the statistical formulation of Barge & Pellat (1990, 1991; see also Hornung et al. 1985). 
Our numerical solutions to the coagulation equation agree with analytical solutions for 
three standard test cases. Our results also agree with WS93's simulation of the formation of 
the Earth at 1 AU. The present models do not incorporate fragmentation of bodies during 
collisions. We will include fragmentation in a separate paper. 

We have considered two simple cases of planetesimal evolution in a 6 AU wide annulus 
centered at 35 AU. Models without velocity evolution invariably produce several large 
bodies that accrete practically all of the material in the annulus. The runaway growth in 
these simulations occurs without dynamical friction or gas drag; it is a direct consequence 
of gravitational focusing. The time required to produce a runaway body in our models 
scales inversely with the initial mass of the annulus and with the initial radii and velocities 
of the planetesimals. For bodies with tq = 80-8000 m and Cinu = 10'^, our simulations 
produce runaway growth in 100 Myr for annular masses of roughly 10-30 Mq. The time 
scale for runaway growth increases to 700-2000 Myr for e = 10~^. A minimum mass solar 
nebula with S oc R~^/'^ contains 7-15 Mq in a 6 AU wide annulus centered at 35 AU. These 
models thus reach runaway growth in a minimum mass nebula on time scales comparable 
to the maximum formation time scale for Neptune (~ 50-100 Myr, p^issauer et al. 1996|; 
Pollack et al. 1996]) . Runaway growth on much shorter time scales, ~ 10 Myr, requires 



annular masses that far exceed the minimum mass solar nebula, ~ 100 in a 6 AU wide 
annulus for Cinit = 10^^ to 10^'^. 

Models with velocity evolution produce runaway growth on a much wider range of time 
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scales compared to constant velocity calculations. First, dynamical friction and viscous 
stirring dominate the evolution of models with relatively large (8 km) initial bodies. The 
velocities of these bodies thus increase as collisions produce more massive objects. This 
velocity increase delays runaway growth by factors of 2 or more compared to constant 
velocity evolution. The delay in the runaway growth time increases with increasing Cinu- 
In contrast, collisional damping dominates the evolution of models with smaller (80-800 
m) initial bodies. These bodies 'cool' until the largest objects have radii of 10-20 km. 
Dynamical friction and viscous stirring then 'heat up' the small bodies but this heating 
is small compared to the velocity increases of the tq = 8 km models. For Cinu = 10~^, 
collisional damping enhances collision rates and decreases the time scale to runaway growth 
by factors of 4-12 compared to constant velocity calculations. Our simulations of minimum 
mass solar nebulae with Cinu = 10~^ reach runaway growth on time scales of 20-40 Myr for 
80 m initial bodies, 50-100 Myr for 800 m bodies, and 75-250 Myr for 8 km bodies. These 
time scales increase by factors of 2-4 for emit = 10~^. 

The formation of runaway bodies in constant velocity simulations is surprising. 
Previous simulations demonstrate a need for dynamical friction, which decreases the 
velocities of the most massive objects ([Wetherill fc Stewart 1989| , 1993; Parge fc Pellat 



imi 1991; lOhtsuki 199^ ; [Ida fc Makino 199^4 1992b; [Kokubo fc Ida 19961 ). Gravitational 



focusing then allows these bodies to sweep up lower mass bodies very rapidly. Kokubo & 
Ida (1996) summarize necessary and sufficient conditions for runaway growth and show that 
the ratio of the maximum mass of planetesimals to their mean mass increases dramatically 
during runaway growth. Our constant velocity simulations satisfy the "necessary" condition 
for runaway growth, 

2Rh 

but do not meet the "sufficient" condition, 

f<0, (3) 

am 

because our velocities are constant with mass, dV/dm = 0. Nevertheless, the maximum 
mass, Mmax, of each constant velocity simulation increases much more rapidly than the 
mean mass, < rrii > (Figure 14). In Figure 14(a), the ratio Mmax/ < f^i > increases slowly 
during the orderly growth phase and then increases rapidly when several isolated bodies 
begin to accrete most of the mass in the annulus. This runaway growth is not as extreme as 
that seen in models with velocity evolution. We derive M^ax/ < rrii > ~ 10^ in constant 
velocity models compared to Mmax/ < rrii > ^ 10*^-10^° in models with velocity evolution 
(Figure 14(b)). Dynamical friction is responsible for the larger increase in Mmax/ < frii > 
in models with velocity evolution. 



- 19 - 



Although the constant velocity simulations are artificial, they are a useful guide for 
planetesimal growth in the outer solar system. Our results show that coUisional damping 
by small bodies can overcome viscous stirring and keep particle velocities low at a = 35 
AU. The relative efficiency of coUisional damping should increase with a, because the 
collision rates decrease more rapidly with a than do the damping rates {tcoii oc a~^/^ vs. 
Tdamp OC a~^/^). The situation in the outer solar system differs markedly from conditions 
at small a. In our simulations at 1 AU, bodies with ^ 80-800 m grow to 10 km objects 
in ~ 1000 yr. As in the 35 AU simulations described above, coUisional damping reduces 
the velocities of the small bodies by a factor of ~ 2 in 1000 yr. In contrast to our 35 AU 
models, viscous stirring and dynamical friction act quickly to increase velocities once larger 
bodies are produced at 1 AU. Aside from the longer time scale to reach runaway growth, 
this evolution then follows closely the simulations in Figures 1-2. 

With these considerations in mind, we suggest a modest modification to the sufficient 
condition for runaway growth, 

f<0, (4) 

am 

instead of the condition in Eqn. (3). This condition maintains the relative growth rate 
necessary for runaway in a three-dimensional system (see Kokubo fc Ida 1996 ) 



1 dm , ^1 /q . , 

— — oc , 5 
M dt ' ^ ^ 

when the gravitational focusing factor is large (e.g., Vij ^ K.ij)- 

In addition to this new criterion for runaway growth, our simulations show that 
Neptune and Pluto can grow in parallel at a ^ 30-35 AU. For a minimum mass solar 
nebula, previous calculations indicate that Neptune reaches its current size in no longer 
than 50-100 Myr ( [ip 1989| ; |Lissauer et al. 1996| ; [Pollack et al. 1996[ ). In our models of 



minimum mass solar nebulae at 35 AU, 800 m objects grow to Pluto-sized planets on similar 
time scales. However, recent observations indicate a much shorter lifetime, ~ 5-10 Myr, 
for gaseous disks surrounding nearby pre-main sequence stars ( [Sargent fc Beckwith 1993| ; 
Strom et al. 1993| ; [Russell et al. 1996| ; see also [Fernandez 1997[) . If our solar system evolved 



on a similar time scale, the formation of the gas-rich outer planets requires solar nebulae 
with masses 2-5 times larger than the minimum mass solar nebula, Mmin ([Lissauer et al. 



[1996[ ; [Pollack et al. 19^B[ ). Neptune then attains its current size in 5-10 Myr. At 35 AU, 



objects reach 1000 km sizes in 10-20 Myr for Mq 2-3 Mmin and in 5-10 Myr for Mq 
5 Mmin assuming the initial bodies have small masses and eccentricities. Fragmentation 
should not change these conclusions unless coUisional erosion can prevent the formation of 
1 km bodies from smaller building blocks. 
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Our Kuiper Belt simulation produce Pluto-sized objects on reasonable time scales for 
other plausible solar nebula models. In particular, Cameron's (1995 and references therein) 
detailed disk models for the protosolar nebula have a shallower density density distribution, 
E = So(a/ao)~^, than the minimum mass model we considered above to derive our success 
criteria (see also [Ruden fc Lin 1986^ [Ruden fc Pollack 1991| ). Cameron's model contains ~ 
100 M0 in our 6 AU annulus. The growth time for 1000+ km objects scales simply from 
results with smaller disk masses. For models with velocity evolution, we estimate r,, ~ 3 
Myr for tq = 80 m, Xr ~ 8 Myr for tq = 800 m, and ~ 15 Myr for tq = 8 km. Although 
all of these models can produce Pluto-sized objects before Neptune reaches its final mass, 
models with ro = 8 km produce more KBOs with ~ 100-300 km as outlined in Sec. 3.3. 

These results - together with recent dynamical calculations of Malhotra (1993, 1995, 
1996) - suggest a self-consistent picture for the formation of Pluto-Charon in the outer 
solar system. In this picture, Pluto and Charon begin as ~ 1 km planetesimals at a ~ 
35-40 AU and grow to their present sizes on a time scale of 10-100 Myr. Both objects are 
runaway bodies more massive than the bulk of the planetesimal mass distribution and have 
low velocities due to collisional damping and dynamical friction. At a somewhat smaller 
semi-major axis, a ~ 25-30 AU, Neptune accretes its current mass in 10-100 Myr and 
migrates radially outward through the protosolar disk during the late stages of giant planet 
formation ( [Malhotra 1993| ; [Pollack et al. 1996| ; [Fernandez fc Ip 1996| ; see also [Ipatov 1989[ , 
1991). During this outward migration, Neptune captures Pluto-Charon in the 3:2 resonance 
( [Malhotra 1995[) . This capture should effectively end further growth of Pluto-Charon, 
because the orbital elements increase to e ~ 0.2 and i ~ 10° on a short time scale, ~ 10 
Myr, inside the resonance ( [Malhotra 1995[ ). High velocity collisions within the resonance 
should also hinder growth as in our large e models. Neptune also captures other KBOs 
at a ^ 35-40 AU into the 3:2 and other resonances ( [Malhotra 199(j[ ; see also [Jewitt et at. 



[1996[) . Further growth of these objects is also slowed due to rapidly increasing velocities 
( [Malhotra 1^ , 1996; [Morbidelh et al. 1995[ ; pavis fc I-'arinella 1997[ ). As long as collisional 
erosion does not decrease significantly the radii of captured KBOs, this sequence of events 
accounts for the general aspects of the mass distribution and orbital elements of observed 
KBOs in a simple way. 

This picture for the formation of Neptune, Pluto-Charon, and KBOs differs from those 
of Stern (1995, 1996) and Stern fc Colwell (1997a,b), who studied KBO evolution at 40-70 
AU. Stern fc Colwell (1997) concluded that the formation of Neptune and KBOs, including 
Pluto-Charon, requires ~ 100-1000 Myr due to the low collision rates at a ^ 30-50 AU. 
Stern fc Colwell's (1997a; also ^tern 19"95[ , 1996) time scale for KBO formation is a factor 
of 10-20 longer than our time scale. Their results also conflict with the more detailed gas 
dynamic calculations of Pollack et al. (1996). Although the exact origin of this discrepancy 
is unclear, we suspect mass spacing, collision rates, and (less probably) fragmentation may 
be responsible. Our mass spacing, 6 = 1.4, should produce more accurate estimates for 
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the growth time than 5 = 2 ( Stern fc Colwell 1997 ) or 5 = 4 (|Stern 1996| ). We estimate 
10%-20% delays in runaway growth for 5 = 2 and expect very long delays for 5 = 4. In 
addition, our collision rates include a factor, l3coih that accounts for the gaussian distribution 
of impact velocities ( [Greenzweig fc Lissauer 1992^ see Eqn. A14). Neglecting this factor 
delays runaway growth by a factor of 3 (see also |WS93| ). Stern (1995) does not include this 
factor in his cross-section (see his Eqn. 5) and thus derives much longer growth times. As 
noted above, our neglect of fragmentation encourages runaway growth. We do not think 
that including this process should delay runaway growth by another factor of 3, for reasons 
outlined earlier, and plan to test this suspicion in our next paper. 

However the theoretical issues may be resolved, several consequences of our accretion 
models can be tested with additional observations. All models that reach runaway growth 
produce 10-100 Pluto-sized objects with radii r^ax ~ 1000-2000 km and a roughly 
power-law mass distribution with a maximum radius at ~ 0.1-0.2rmax (see Figures 8-10). 
If Pluto and Charon are runaway bodies produced at a 35-40 AU, there should be 
several additional "Plutos" with similar orbital elements. Malhotra (1995) reached a similar 
conclusion and noted that recent searches have not excluded the possibility of several 
additional Plutos in 3:2 orbital resonance with Neptune. The accretion models also predict 
a cumulative power law distribution, Nc oc r~^, for objects with r ^ 100-200 km. By 
analogy with WS93, this shape is fairly independent of fragmentation as long as collisions 
produce overall growth instead of disruption. The best fit to the observed distribution is 
shallower than expected, Nc oc r~^, but the data are not yet accurate enough to preclude 
our model prediction (|Jewitt et al. 1996|) . Larger surveys will provide a better test of this 
prediction ( [Luu et al. 1997| ). 



Finally, our results suggest that KBO formation is likely in other solar systems. KBOs 
can grow in the dusty disks that surround many nearby main-sequence stars, if the disk 
masses are within an order of magnitude of the "maximum" disk masses for a Lyr, a PsA, 
and f5 Pic ( [Backman fc Paresce 199B ). Stars with disk masses near the minimum dust 
masses of these well-studied A-type stars probably have few, if any, KBOs, but could have 
many objects with ~ 1-10 km. Observations of nearby pre-main sequence stars also 
indicate substantial masses, ^ 10^ circumstellar disks with radii of 100-1000 

AU (e.g., ^argent fc Beckwith 1993| ; see also [Beckwith fc Sargent 1996|) . These data imply 
Mo ~ 10-100 M0 in the Kuiper Belt. With formation time scales of 1-10 Myr, KBOs can 
grow in massive circumstellar disks during the pre-main sequence contraction phase of a low 
mass star (see Kenyon fc Hartmann 1995| and references therein). In a less massive disk, 
KBOs grow while the central star is on the main sequence. We expect no KBO formation in 
circumstellar disks with very low masses, ^ 10 M^. These disks can produce 1-10 km 
objects at a 30 AU, unless fragmentation prevents growth of icy bodies. We will explore 
this possibility in a subsequent paper. 
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A. APPENDIX 



A.l. Overview 

Our evolution model follows procedures developed for other planet formation 
calculations, including Safronov (1969), Greenberg et al. (1978, 1984) and WS93. We 
assume planetesimals are a statistical ensemble of masses with a distribution of horizontal 
and vertical velocities about a single Keplerian orbit. We consider a cylindrical annulus of 
width, Aa, and height, if, centered at a radius a from the Sun. Particles in the annulus 
have a horizontal velocity, /ii(t), and vertical velocity, fj(t), relative to an orbit with mean 
Keplerian velocity, Vk (see [Lissauer fc Stewart 19931) . These velocities are related to the 
eccentricity, Cj, and inclination, ij, through 

= + lsin\)Vl (Al) 

with 



K 



(A2) 



and 



1 



sin ijVi 



K 



(A3) 



We approximate the continuous distribution of particle masses with discrete batches having 
an integral number of particles, riiit), and total masses, Mj(t) ( [WS93| ). The average mass 
of a batch, mi{t) = Mi{t)/ni{t), evolves with time as collisions add and remove bodies 
from the batch. This procedure naturally conserves mass and allows a coarser grid than 
simulations with fixed mass bins ( [Wetherill 1990| ; see also [Ohtsuki fc Nakagawa 1988| ; 
Kolvoord fc Greenberg 1992| ). 
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To follow the evolution of particle number, we solve the coagulation equations for all 
mass bins, k, during a time step, 6t, 

Suk = St [eijAijUiUj - rikAikUi] - Suk^gd (A4) 

SMk = St [eijAijUinjUik - nkAikUimk] - rUkSuk^gd (A5) 

where Aij is the cross-section, eij = 1/2 for i = j and eij = 1 for i ^ j. The three terms in 
A4-A5 represent (i) mergers of rrii and rrij into a body of mass = mj + rrij, (ii) loss of 

through mergers with other bodies, and (iii) loss of by gas drag. This treatment 
assumes (i) that each body can collide with every other body and (ii) that bodies do not 
fragment during collisions. Assumption (i) is correct for all but the very largest bodies, 
which become isolated from one another as their orbits circularize due to dynamical friction 
(see below). We correct equations A4-A5 for this effect by calculating the "gravitational 
range" of the largest bodies - Rg^i = KiaRn^u^i^ + 2aei ( |WS93|) - where Ki = 2v^ and 



RH,ij = [{'iTT'i + ruj) /^MqY^^ is the mutual Hill radius. As in WS93, the isolated bodies 
are the largest bodies that satisfy the summation, Sl^"^ ^iRg,! > ^o-- Assumption (ii) 
is rarely correct, because all collisions produce some debris unless the relative velocity of 
the two particles is very low (see, for example, [WS93|) . In this paper, we concentrate on 



planetesimal growth and assume that all collisions result in mergers. We will consider the 
effects of fragmentation in a separate paper. 

To calculate the appropriate index k for a specific collision between batches i and j, 
we first calculate a fixed grid of masses, m;, for Z = 1 to Nmax and S = mi+i/mi. The mass 
spacing, S, is constant throughout a calculation; N^ax increases with time as more batches 
fill with particles. When a collision produces bodies with m^, we augment either batch I 
when rrik < y/minii^i or batch / + 1 when m^, > ^m^m^+i. A complete cycle through all 
mass batches produces new values for and M^, which yields new values for the average 
mass per bin, ruk = Mk/uk- This process conserves mass and provides a good description 
of coagulation when S is small (see below) . 

Besides collisions, several processes contribute to the velocity evolution of growing 
planetesimals, including dynamical friction, gas drag, and viscous stirring. We assume that 
all collisions between mass batches conserve the horizontal and vertical components of 
kinetic energy, E^^i = mihf/2 and E^^i = mivf/2. The change in the two components of 
kinetic energy due to collisions is 

= -1/2 [Sukimkhl) = Sn,{m,h1) + Snj{mjh])\ (A6) 
SEl\ = -1/2 ISnkinikvl) = Sni{miV^) + Snj{mjV^)] (A7) 
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for each pair of collisions between and rrij. In these expressions, Sui < represents 
the change in rii due to collisions with particles in batch j. Batch k loses kinetic energy 
due to collisions with other batches (e.g., 6nk > 0). We also calculate the evolution of hi 
and Vi due to gas drag ( |Adachi et al. 1976D and collective interactions, such as dynamical 



friction and viscous stirring, using a statistical treatment of the appropriate Boltzmann 



and Fokker-Planck equations (|Hornung et al. 1985| ; see section A. 3 below). The complete 



change in the horizontal and vertical kinetic energies is thus: 



6E,,, = 6E:% + 5El% + + 5El, (A9) 

where the superscripts "gd" (gas drag), "in" (inelastic), and "Ir" (long range elastic 
collisions, such as viscous stirring and dynamical friction) refer to a specific type of velocity 
evolution outlined below (see also [Barge &: Pellat 19901 , 1991; [WS93|) . 



We solve the complete set of evolution equations, A4, A5, A8, and A9, using an explicit 
method that automatically prevents large changes (> 0.1%) in the dynamical variables - n^, 
Mj, hi, and Vi - by limiting the time step. As in WS93, we require integer values for rij and 
5ni. Section A. 4 compares our numerical procedures with analytic results from Wetherill 
(1990; see also Phtsuki fc Nakagawa 198^ ; [Ohtsuki et al. 1990|) . Section 2 of the main 



text compares calculations at 1 AU with results from WS93. In both cases, our procedures 
reproduce the expected results. Before describing the analytic results, we first describe in 
detail our treatment of the collision rates (section A. 2) and the velocity evolution (section 
A.3). 



A.2. The Collision Rate 

Approximations to the collision rates between planetesimals are in the spirit of kinetic 
theory, where the number of collisions is the product of the local density, the relative 
velocity, and a cross- sect ion. WS93 express the number of collisions between a single body, 
mj, and all of the bodies, mj, as 

nc,^J = a^u {jjf^^) V^J F9,ij i^i + rjf 5t (AlO) 

where and rj are the radii of the two bodies; V^j = + is the relative velocity; 
and Fg^ij is the gravitational focusing factor. The constant factor acou accounts for the 
gaussian distribution of particle velocities and the difference between the collision frequency 
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of particles on Keplerian orbits and those in a box ( Greenzweig fc Lissauer 1992| ; |WS93|) . 
The relative velocities and scale height depend on the individual particle velocities ( [WS93|) : 



(All) 



V,, = (h^ + vf + e + v')'/' 



(A12) 



where Q is the Keplerian angular frequency. The total number of collisions for rrii is simply 
nific] the cross-section appropriate for equations A4-A5 is then: 



A. 



OicoU 



1 



4 H a Aa 



(A13) 



We consider two approaches to compute the gravitational focusing factor, Fg^ij. In the 
first case, we follow WS93 and set 



WS,ii 



Eij {1 + /3, 



coll ' 



(AM) 



where V^^j = 2G{mi + mj) / (r^ + rj) is the mutual escape velocity. The extra factors account 
for the gaussian distribution of impact velocities {jScoii, [Greenzweig fc Lissauer 1992|) and 
the deviations from two-body focusing at low relative velocities {Eij, preenzweig fc Lissauer 
|1990|) . We adopt WS93's prescription for the variation of Pcoii as a function of the relative 
velocity in Hill units, VH,ij = Vij/{RH,ijVK/a): 



coll 



2.7 

1.0 + 1.7(\/h,^,-1) 
1.0 



1 < Vma. < 2 



V, 



H,ij 



H,ij 
< 1 



(A15) 



and we set 



E., 



(A16) 



where Ek = Jq^'^ a/I — fc^sin^^ d9 for k'^ = 3/[4(l -|- sin i/e)] ( Greenzweig fc Lissauer 1990 ; 
Greenberg et al. 1991| ). 

At very low velocities {yH,ij < 2.3), we adopt the two-body collisional cross-sections of 
Greenberg et al. (1991): 
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2B,ij 



i 0.5(1 +Ky^,r^(f) (ft) 



VH,ij < 2.3 



VH,ij < 2.3,VH,ij < Vh, 



(A17) 



crit 



where Rt = a[{mi + rrij)/! M©]^/^ is the Tisserand radius, Vr = LlilAar is the Tisserand 
velocity, Aa^ = 2.5Rh is the half width of the feeding zone, VH,ij = Vij/Vn is the relative 
vertical velocity in Hill units, and VH,crit < 0.7 sin (0.9 [(rj + rj)/ Rn^ijY^'^)- 

In our second approach to gravitational focusing, we modify the piecewise analytical 
approximation of Spaute et al. (1991): 



l+/3eo.(^) 

42.4042 (^) 
10706.916 



Vij > 0.032K,ii 

0.01 < Vij/Ve,ij < 0.032 

V^j/Ve,^, < 0.01 



(A18) 



with Pcoii as defined above. These expressions are continuous and serve as a check on the 
more detailed expressions for Fws,ij- For very low velocity encounters, we use the two body 
cross-sections for F2B defined above. 

In these approximations to the cross-section, the transition from Fws or Fs to F2B at 
VH,ij ~ 2.3 is not smooth. To affect a smooth transition, we set 



[1 — X2B,ij) Fws,ij + X2B,ij F2B,ij 
[1 — X2B,ij) Fs^ij + X2B,ij F2B,ij 



(A19) 



where 



X2B,' 





Q.^{VH,ij - 2.3) 
1 



VH,ij > 3.3 

1.3 < VH,ij < 3.3 

VH,ij < 1.3 



(A20) 



A. 3. Velocity Evolution 



As noted above, kinetic models approximate planetesimal orbital elements as a mean 
square random velocity, V^^ ( ^afronov 1969^ [Stewart fc Wetherill 1988| ; [WS93| and references 
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therein). We divide this velocity into horizontal, hi, and vertical, fj, components that 
are related to Vi, e, and i (see equations A1-A3). Hornung et al. (1985) derive analytic 
expressions for the time evolution of planetesimal velocities using a kinetic approximation 
to average over the velocity distribution function. WS93 reformulate some of these results 
in terms of the eccentricity and inclination, which we adopt here for simplicity (see also 
Stewart fc Wetherill 1988| ). We calculate velocity changes due to (i) gas drag, which 



decreases particle velocities and causes particles to spiral in through the disk; (ii) dynamical 
friction from elastic collisions, which transfers kinetic energy from larger to smaller bodies; 
(iii) viscous stirring from elastic collisions, which taps the solar gravitational field to increase 
the velocities of all bodies; and (iv) collisional damping from inelastic collisions (see also 
Hornung et al. 1985| ; [Barge fc Pellat 199q , 1991; [Ohtsuki 19921) . 



The time evolution of the eccentricity and inclination for long-range, elastic encounters 
is ([WS93l Appendix C) 

^ = E ^ K + el {Jr + AJe) (A21) 

, •2 j=N ^ 

if =E^K + -.)^?^. (A22) 



for viscous stirring and 



J 2 j=N ^ 

if = E ^ (^.S' - m,eD {Kr + 4Ke) (A23) 

, •2 j=N ^ 

¥=T.^ - (A24) 



dt 



for dynamical friction. In these expressions, Cj, e^, Zj, and ij are the eccentricity and 
inclination of each body, = {if + «^)/(e^ + e^) is the ratio of inclination to eccentricity, 
and Clr = 16G^pj(lnA + 0.55)/V^(e^ + e|)'^/^ is a function of the density of particles in 
batch j and the relative horizontal velocity of the mass batches ([WS93| ). The functions 
Jj., Je, Jz, Kr, Kg, and K.^ are definite integrals that are functions only of Pij ( [Hornung ei| 
\al. 1985[ ; Parge fc Pellat 1990[ , 1991; [WS93[ ; phtsuki 1992[ describes a similar approach to 



velocity evolution). 

The time evolution of e and i due to collisional damping is 

j=i (-i 

if = E ^ (^.4 - m,e^ - (m, + m,)eD (/. + 4J,) (A25) 
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"-''''171.. 7, \ ^ 



m,i '-'m / -2 -2 ( , \ -2 



rrii + m,)%\) h (A26) 



where Cin = acoiiUjPjVijFg^ij{ri + rjY and pj is the mass density ( pornung et ai 1985| ; 



see also [Ohtsuki 1992|) . We include terms from the collision rate, acoiu ^ij, and Fg^ij, for 



consistency. The integrals, Irif^ij) hif^ij), and Izif^ij), are listed in the Appendices of 
Hornung et al. (1985). We integrate these expressions numerically. The integral often 
diverges; we set Iz = I — + h) to avoid these divergences. 

In addition to dynamical friction and viscous stirring, we also consider velocity 
evolution due to gas drag. Gas drag reduces the velocities of all mass batches and also 
removes material from each mass match. The inward drift of material is ( [Adachi et al. 



— = 2(0.97e + 0.64Z + r//V^) - , (A27) 



where rj is the gas velocity relative to the local Keplerian velocity, Vk- We adopt = 60 m 
s~^ for calculations at 1 AU ([WS93| ) and 77 = 30 m s~^ for calculations at 35 AU ([Adachi e\ 



ai. 19761). The characteristic drift time is 




-0 = ^ T7^ ^ \Tk. (A28) 



where Cd = 0.5 is the drag coefficient, pg = 1.18 x 10 ^ (a/1 AU) ^^^'^ is the gas density 
(Nakagawa et al. 1983), and is the orbital period (see [Adachi et al. 197"^ , [WS93[ ). To 



simulate the disappearance of gas in the protosolar nebula, we decrease the gas density with 
time: 



Pg{t) = Pgfl e 



(A29) 



with 



1.18 X 10^9 g cm-3 (a/1 K\])-^^/^ (Mo/Mmin) ( [Kuden fc Lin 1986[ ; |Ruden fc 



Pollack 1991[ ; [WS93| ). The radial decrease of the gas density follows models for minimum 
mass solar nebulae; the mass dependence allows the density to scale with the mass of the 
annulus. 

The number of bodies lost from the calculation at each time step depends on their 
effectiveness at crossing Aa. We set the number of bodies lost from a batch as: 



5n, 



/Aa\ (5V 

Hi \ a J \To^ 



gd,i 



(A30) 
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This expression is used in equations A4-A5. 

Finally, we adopt Wetherill & Stewart's (1989) expression for velocity damping due to 
gas drag: 



dV, -ttCt 



dt 



2m,- 



(A31) 



where Co = 0.5 is the drag coefficient and Vg = {Vi{Vi + r]))^^"^ is the mean relative velocity 
of the gas. 

We convert the differential equations, A21-A26 and A31, into a kinetic energy form in 
two steps. We use Pi to derive the appropriate horizontal and vertical components of the 
velocity, V^, in Eqn A31. Equations A2-A3 similarly yield 5hi and 5vi in terms of 5ei and 
5ii, where 



and 



These substitutions yield: 



5ef = St 



5if = 5t 



d^ 
dt 

df 
dt 



TrCppgVg^rf ^ 
2mi{l + O.SPi) 



(A32) 



(A33) 



(A34) 



{Svfr = 0.8 Sh' 



(A35) 



{SKrr = -Vl 6t 



(Svff = - sin^i 5t 



de^ ■ 
dt 



dt 



lr\2 



^Vi5t 



de: 



dp"^ 



{5vtf = - sin^i Vl St 



dt 

'di^ 



di^ 1 



dt dt 

We multiply these relations by rrii for substitution into equations A8-A9. 



(A36) 
(A37) 
(A38) 
(A39) 
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A. 4. Tests of the evolution code 

To test the validity of our numerical techniques, we compare our results with several test 



cases see 



Uhtsuki fc Nakagawa 198^ ; [Ohtsuki et al. IQQQ ; [Wetherill 19"90|) . The coagulation 



equations, A4-A5, have analytic solutions for three simple forms of the cross-section, 
Aij. Smoluchowski (1916) first solved the coagulation equation for Aij = ac = constant. 
Trubnikov (1971) described solutions for A^j = /?c(^i + ttlj) and Aij = '-jcfriimj. Wetherill 
(1990) identified an inconsistency in Trubnikov's results for Aij = '-ycfni^j and showed that 
this cross-section produces runaway growth. Tanaka & Nakazawa (1994) verified Wetherill's 
new solution and placed limits on the validity of the coagulation equation during runaway 
growth. 

The analytic solutions to the coagulation equation provide rigorous tests of numerical 
methods. Two simple cases, Aij = ac and Aij = (3c{Tni + rrij), do not lead to runaway 
growth, but they test the ability of numerical codes to reach a target mass at a specified 
time. They also yield estimates for mass conservation over many time steps. Numerical 
solutions for Aij = jc'^imj are challenging, because runaway growth requires a careful, 
automatic procedure for changing the time step. In all three cases, the time lag between the 
analytic and numerical solutions depends on the mass ratio between consecutive batches, S 
( [Wetherill 1990|) . These tests thus yield a quantitative measure of the largest allowed value 
for 6 (lOhtsuki fc Nakagawa 1988| ; [Ohtsuki et al. 1990| ; [Wetherill 1990|) . 



To compare our numerical results with analytic solutions, we follow conventions 
established by Ohtsuki et al. (1988) and Wetherill (1990). For Aij = and 
Aij = Pcirrii + rrij), we plot log Nimf as a function of log m^. We evaluate log Nk as a 
function of log and the fractional mass in the small body swarm and the runaway body 
as a function of a dimensionless time for Aij = '-jcTrtinij. The numerical calculations do not 
have regular mass intervals, so we calculate N = 6N/6m, where 

6N = Nk + {Nk+i + Nu-i)/2 (A40) 



3m = rrik+i - rrik-i. (A41) 

Wetherill (1990) describes analytic solutions for each cross-section in detail. The 
Aij = ac case has the simplest solution. If uq is the initial number of particles with mass 
rrii, the number of bodies with mass = knii at a time, t, is 

rik = riofil - ff-' (A42) 
where / = 1/(1 + ?7i/2) and rji = acUot is the dimensionless time (see also pilk fc Takahashi] 



1979| ; phtsuki fc Nakagawa 19881 ). The solution for Aij = (3c{mi + rrij) has a similar form: 
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nu = ^0^/(1 - ff^'e^^^'-f^ (A43) 

where / = e~^'^ and 772 = (3cnot. In both of these expressions, / is the fraction of bodies 
with irii that have yet to undergo a colhsion at time t. 

Models with Aij = '-yjriimj lead to runaway growth when rj-^ = 'ynot = 1 (IWetheriH 



19901 ; Parge fc Pellat 1990^ P?anaka fc Nakazawa 19941) . The number distribution for r/3 < 1 



IS 

This solution fails to conserve mass for 773 > 1; a single runaway body then contains most 
of the total mass. The mass of the runaway body for 173 > 1 is ([Wetherill 1990| ; P?anaka S 
Nakazawa 1993| , 1994): 



mn = noe 



Figure 15 compares our results for 6 = 1.25 with the analytic solution for Aij = ac- 
The agreement is good and again improves as 5 decreases. Our results for 5 = 1.4-1.6 
are consistent with the analytic solution, although we have too few mass batches to make 
reliable comparisons when 771 < 10 — 20. We did not attempt numerical models for 6 = 
1.6-2 (the maximum allowed), but we expect that these will produce satisfactory results for 
7] > 100. 

Figure 16 shows results for A^j = Pd'^i + f^^j) and 6 = 1.25. The agreement between 
our calculation and the analytic solution is quite good and improves as 6 increases. We 
find a slight excess of low mass bodies in our numerical results compared to the analytic 
solution. Wetherill's Figure 4 contains a similar excess. The peak of our normalized number 
distribution lags the analytic result by 1.4%. This lag decreases with 6 and is < 1% for 6 = 
1.10. 

Figures 17 and 18 summarize our results for A^j = 'yciTiirrij and 6 = 1.25 for 77 < 1. The 
numerical solution follows the analytic model very closely for 7/3 < 0.95 and then begins 
to diverge at large masses as rj approaches unity (Figure 17). The numerical model begins 
runaway growth at 773 = 1.012 and lags the analytic model by 1.2%. The numerical runaway 
begins much closer to the predicted result, 773 = 1.005, for 6 = 1.08. Larger values for 6 
produce runaways that are delayed by much longer factors. The lag is 2.7% for 6 = 1.4 and 
8.7% for 6 = 2. Wetherill (1990) quotes similar results for his numerical models with 6 = 
1.07 and 6 = 1.25. 

Figure 18 describes the evolution of the runaway growth model for 6 = 1.08 and 773 > 
1. The calculated mass distribution initially lags the analytic result by less than 1% for 773 
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marginally larger than 1 (see also |WS93 ) but matches the analytic result almost exactly at 
rjs = 1.05 (Figure 18; left panel). The calculation continues to match the analytic result 
until rjs ^ 5. The right panel of Figure 18 plots the mass of the runaway body for 773 > 1. 
The calculated mass agrees with the analytic prediction, Eqn. A45, to 1% or better for all 
T]3 > 1. Models with S = 1.25 have greater difficulty reaching large 773 due to their poorer 
mass resolution. These models have larger at high masses, which reduces the time step 
considerably compared to models with small 6. The calculation then requires a significant 
amount of computer time and does not agree as well with the analytic predictions. Our 
models with 6 > lA fail to reach rj^ ^ 1.1 if we maintain our criterion of small Srif^ per 
time step. Relaxing this criterion allows reasonable time steps but produces very poor 
agreement, ^ 20%, with the analytic solution. 

These results confirm our limits on Srik for the Kuiper Belt simulation described in 
the main text, ^ 0.1% per time step, for 6 = 1.08-1.4. Models with larger S fail to follow 
growth properly unless the time steps are unreasonably small. 
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Table 1. Basic Model Parameters 



Paraiiietor 


Syml)ol 


1 Ai; Models 


35 Ai; Models 


Width of annulus 


Sa 


0.17 AU 


6 AU 


Initial Velocity 


Vo 


4.7 m s~^ 


4.5-45 m s~-^ 


Particle mass density 


Po 


3 g cm~^ 


1.5 g cm~^ 


Relative gas velocity 


V 


60 m s~^ 


30 m s-i 


Time Step 


St 


0.5 yr 


5-250 yr 


Number of mass bins 


N 


100-150 


64-128 


Mass S])aciiig of l)ins 


(S 


< 1.20 


1.10 



Table 2. Model Results at 1 AU^") 







5 = 1.25 






5 = 1.40 




Time (yr) 


Tmax (km) 


m{r,„ax) (kg) 


ij'rriax) 


Tmax (km) 


m{r,„as) (kg) 


{j max ) 


5.0 xlO^ 


19.8 


9.7 xlO^^ 


3 


20.9 


1.2 xlO^o 


1 


1.0 xlO^ 


513.3 


1.7 xl024 


3 


492.3 


1.5 Xl024 


1 


2.5 xlO^ 


1167.3 


2.0 xlO^^ 


1 


1203.0 


2.2 xlO^s 


1 


5.0 xlO^ 


1540.8 


4.6 xlO^s 


3 


1515.7 


4.4 xlO^s 


1 


1.0 xlO^ 


1890.8 


8.5 xlO^s 


2 


1746.7 


6.7 xlO^s 


3 


1.5 xlQ-'^ 


1918.3 


9.3 xlO^^'^ 


5 


2382.7 


1.7 xlO^fi 


1 



•^"^ These results are for the WS93 prescription of gravitational focusing, equations 
A14-A16, as summarized in the main text. 
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Table 3. Model Results at 1 AU^") 







5 = 1.25 






5 = 1.40 




Time (yr) 


Tmax (km) 


m{r„,ax) (kg) 


ij'rriax) 


Tmax (km) 


m{r,„ax) (kg) 


{j max ) 


5.0 xlO^ 


19.8 


9.7 xlO^^ 


3 


20.9 


1.2 xlO^o 


1 


1.0 xlO^ 


550.8 


2.1 X1024 


2 


492.3 


1.5 Xl024 


1 


2.5 xlO^ 


1167.5 


2.0 xlO^^ 


1 


1223.0 


2.3 xlO^s 


1 


5.0 xlO^ 


1458.6 


3.9 xlO^s 


4 


1735.0 


6.6 xlO^s 


2 


1.0 xlO^ 


1969.0 


9.6 xlO^s 


1 


2174.8 


1.3 xlO^e 


1 


1.5 xlQ-'^ 


2121.0 


1.2 xlO^fi 


1 


2335.0 


1.7 xlO^fi 


1 



•^"^These results are for an adaptation of the Spaute et al. (1991) prescription of 
gravitational focusing, equation A18, as summarized in the main text. 
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Table 4. Model Results at 100 Myr (No Fragmentation, No Velocity Evolution) 



Mo (Me) e ro (m) Nq rg^x (km) rg (km) rmax (km) N{rmax) Tr (MYr) 



1 


10-3 


80 


1.87 xlOis 


0.5 


0.6 


0.9 


1 


2698 


3 






5.61 xlO^"" 


1.1 


1.5 


2.2 


1 


902 


10 






1.87 xlO ° 


3.9 


5.4 


6.7 


1925 


270 


30 






5.61 xlO^*^ 


NA 


NA 


8500 


1 


90 


100 






1.87 xlO" 


NA 


NA 


18500 


1 


27 


1 




800 


1.87 xlOi2 


1.4 


2.2 


2.7 


1 


2340 


3 






5.61 xl0^2 


2.2 


3.0 


4.2 


4 


780 


10 






1.87 xlO^-' 


5.5 


8.5 


10.7 


341 


234 


30 






5.61 xlO^^ 


NA 


NA 


18000 


2 


78 


100 






1.87 xlO^^ 


NA 


NA 


16700 


1 


23.5 


1 




8000 


1.87 xlO^ 


13.3 


16.9 


20.7 


140 


753 


3 






5.61 xlO^ 


19.2 


26.8 


36.7 


91 


250 


10 






1.87 xlO^° 


NA 


NA 


7585 


5 


75 


30 






5.61 xlO^*^ 


NA 


NA 


9800 


1 


25 


100 






1.87 xlO^^ 


NA 


NA 


8700 


2 


7.5 


1 


10-2 


80 


1.87 xlO^s 


0.5 


0.6 


0.9 


1 




3 






5.61 xlOi'^ 


1.1 


1.5 


1.9 


634 




10 






1.87 xlOi6 


3.1 


4.8 


6.0 


28 


2552 


30 






5.61 xlOis 


9.6 


13.2 


16.4 


323 


850 


100 






1.87 xlO^^ 


33.6 


46.5 


63.5 


5 


255 


1 




800 


1.87 xlOi2 


1.4 


2.2 


2.4 


3195 




3 






5.61 xlOi2 


2.2 


3.1 


3.8 


75 




10 






1.87 xl0i3 


4.4 


6.8 


7.4 


2018 


2521 


30 






5.61 xl0i3 


10.8 


14.8 


18.4 


387 


840 


100 






1.87 xlOi4 


37.6 


51.8 


65.0 


70 


251 


1 




8000 


1.87 xlO^ 


10.1 


11.6 


16.1 


1 




3 






5.61 xlO^ 


11.7 


14.9 


18.5 


8 




10 






1.87 xlOi° 


13.4 


18.5 


23.3 


298 


2146 


30 






5.61 xlOi° 


21.2 


29.3 


39.9 


1 


716 


100 






1.87 xlO" 


59.2 


81.9 


104.2 


18 


215 



Table 5. Model Results (Velocity Evolution, No Fragmentation) 



Mo (Me) 


e 


ro (m) 


No 


r95% (km) 


r-5 (km) 


rmax (km) 


N {Vyriax) 


Tr (MYr) 


1 


10^2 


80 


1.87 X 10^^ 


0.5 


0.8 


0.9 


6 




2 






3.74 X 10^^ 


2.6 


30.4 


126.1 


2 




3 






5 61 X in^^ 


2.7 


33.8 


603.4 


1 


184 


6 






1.12 X 10^^ 


2.9 


44.1 


1072.6 


1 


67 


10 






1.87 X 10^^ 


2.6 


54.1 


1000.5 


1 


32 


17 






3 18 X m^^ 


2.4 


55.7 


1008.5 


1 


18 


30 






5.61 X IQis 


2.7 


74.4 


1059.2 


1 


10 


1 




800 


1.87 X 10^^ 


1.4 


2.3 


2.7 


9 




2 






3.74 X 10^^ 


2.0 


3.2 


3.7 


24 




3 






5 61 X in^^ 


2.8 


5.4 


6.8 


1 




6 






1.12 X 10^^ 


11.2 


47.3 


456.9 


1 


155 


10 






1.87 X 10^^ 


10.6 


51.4 


1007.8 


1 


83 


17 






3 18 X m^^ 


10.6 


68.2 


1022.7 


1 


46 


30 






5.61 X IQi^ 


10.6 


90.7 


1000.3 


2 


25 


1 




8000 


1.87 X 10^ 


11.7 


17.9 


23.9 


8 




2 






3.74 X 10^ 


15.3 


24.1 


32.7 


7 




3 






5.61 X 10^ 


17.3 


30.6 


46.0 


4 




6 






1.12 X 10^° 


27.3 


57.1 


123.0 


1 




10 






1.87 X 10^° 


43.4 


90.6 


462.4 


1 


132 


17 






3.18 X 101° 


48.8 


101.2 


1050.0 


1 


76 


30 






5.61 X 101° 


51.7 


114.2 


1017.7 


2 


44 


1 


10~2 


80 


1.87 X 10^^ 


0.4 


0.6 


0.8 


762 




2 






3.74 X 10^^ 


0.7 


1.2 


1.4 


39 




3 






5 61 X Ifii^ 


1.0 


1.8 


1.9 


1103 




6 






1.12 X 10^^ 


5.8 


46.6 


294.7 


1 


126 


10 






1.87 X 10^^ 


6.3 


60.4 


1004.1 


1 


76 


17 






3 18 X Ifii^ 


6.4 


70.0 


1002.2 


2 


43 


30 






5.61 X 10^^ 


6.4 


82.2 


1070.2 


1 


24 


1 




800 


1.87 X 10^^ 


1.4 


2.2 


2.7 


2 




2 






3.74 X 10^^ 


1.8 


2.8 


3.4 


2 




3 






5 61 X 10^2 


2.2 


3.3 


3.8 


48 




6 






1.12 X 10^^ 


3.1 


4.9 


6.0 


1 




10 






1.87 X 10^^ 


4.4 


7.3 


8.4 


78 




17 






3 18 X m^^ 


13.9 


43.0 


74.1 


4 


130 


30 






5.61 X 10^^ 


23.4 


90.8 


1000.5 


2 


73 


1 




8000 


1.87 X 10^ 


10.1 


13.0 


16.3 


6 




2 






3.74 X 10^ 


10.1 


14.3 


16.4 


303 




3 






5.61 X 10^ 


11.7 


15.3 


18.3 


11 




6 






1.12 X 101° 


13.2 


17.8 


20.7 


68 




10 






1.87 X 10i° 


13.4 


20.8 


26.1 


2 




17 






3.18 X 10i° 


17.0 


25.9 


31.0 


46 




30 






5.61 X 10i° 


28.3 


45.8 


62.3 


6 


167 
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FIGURE CAPTIONS 



Fig. 1. — Results at 1 AU for Mq — 0.667 and S — 1.25. (a) Cumulative mass distribution 
at selected times. The "runaway plateau" forms at ~ 2 x 10"^ yr; it includes 17% of the total 
mass at 5 X 10"^ yr, 18% of the total mass at 10^ yr, and 22% at 1.5 x 10^ yr. (b) Horizontal 
velocity distribution. Viscous stirring increases all velocities with time; dynamical friction 
breaks the runaway bodies and increases velocities of the lowest mass bodies. 

Fig. 2. — Results at 1 AU for 5 — lA. The runaway plateau is more poorly resolved as the 
mass spacing increases, (a) Cumulative mass distribution at selected times. The "runaway 
plateau" includes 12% of the total mass at 5 x 10"^ yr; this component comprises 18% of 
the total mass at 10^ yr and 22% at 1.5 x 10^ yr. (b) Horizontal velocity distribution. 
The horizontal velocities of low mass objects are 5% larger with 6 — 1.4 than with S — 
1.25. Higher mass bodies have velocities 2-3 larger than comparable masses in S = 1.25 
simulations. 

Fig. 3. — Cumulative size distribution for Mq = 10 M^, tq = 8 km, and e = 10~^. The 
eccentricity of this model is constant in time. CoUisional growth is quasi-linear for 50 Myr 
until the largest bodies have Vmax — 50 km. Runaway growth begins when r^ax ^ 100 km; 
these bodies then grow to sizes of 10^ km to 10'^ km in 20-30 Myr. 

Fig. 4. — Radius evolution as a function of initial mass for constant velocity models, with 
(a) To = 8 km, (b) Tq = 800 m, and (c) Tq = 80 m. The eccentricity is e = 10~^. The time 
to runaway growth scales inversely with initial mass, ~ tq (Mo/1M0)-^; To ^ 753 Myr 
for ro = 8 km, tq ^ 2340 Myr for tq = 800 m, and Tq ^ 2700 Myr for tq = 80 m. 

Fig. 5. — Same as Fig. 4, but with eccentricity e = 10~^. (a) ro — 8 km, (b) Tq = 800 m, 
and (c) ro = 80 m. 
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Fig. 6. — Cumulative size distribution for a constant velocity model with 5 = 1.1 for Mq — 
10 M0, ro = 8 km, and e = 10~^. This model produces a cumulative size distribution with 
two distinct power laws: Nc oc r~^'^ for 8 km ^ ^ 300 km and Nc oc r~^^'^ for 300 km 
^ n ^ 5000 km. 

Fig. 7. — Time scales for viscous stirring (solid line) and coUisional damping (dashed line) 
as a function of mass, (a) The two curves show the time scales for interactions only between 
particles of the same mass for a reahstic cumulative mass distribution. The two plots curve 
up at high masses due to the small number of particles. CoUisional damping balances the 
velocity increase due to viscous stirring at ~ 10^^ kg. (b) Time scales integrated over all 
particles during the initial stages, ~ 10 Myr, of a model with velocity evolution. CoUisional 
damping overcomes velocity increases from viscous stirring only for <^ 10^^ kg. (c) As in 
(b) but for a late stage in the growth process, 15-16 Myr. The integrated effects of viscous 
stirring now increase the velocities of all particles. 

Fig. 8. — Model with Mq — 10 M^, Vq — S km, and e = 10~^, with velocity evolution: (a) 
cumulative size distribution, and (b) horizontal velocity as a function of time. CoUisional 
growth is quasi-linear until the largest bodies have r„iax = 50-100 km at 50-60 Myr. Runaway 
growth begins when rmax ^ 500 km at ~ 100 Myr; these bodies then grow to sizes of 10^ 
km in another 50-80 Myr. 

Fig. 9.— Model with Mq = 10 Mq, ro = 800 m, and e = 10~^: (a) cumulative size 
distribution, and (b) horizontal velocity as a function of time. CoUisional growth is quasi- 
linear for 45-50 Myr until the largest bodies have r„iax = 50-100 km. The transition to 
runaway growth requires ~ 10 Myr, when Charon-sized objects form. These bodies grow to 
sizes of 1000 km in another 20 Myr. The velocities of the smallest objects increase with due 
due to viscous stirring. Dynamical friction reduces the velocities of the largest objects. The 
velocity minimum at 3-5 km indicates the batches that contain the largest fraction of the 



-44- 



total mass. 

Fig. 10.— Model with Mq = 10 M0, tq = 80 m, and e = 10"^: (a) cumulative size 
distribution, and (b) horizontal velocity as a function of time. The largest bodies reach rmax 
= 50-100 km in 20 Myr and r^nax = 500 km in only 25 Myr. Runaway growth begins at 25 
Myr and the largest bodies achieve Vmax ~ 1000 km after another 8 Myr. As in Figure 9, 
dynamical friction and viscous stirring increase the velocities of the smallest objects at the 
expense of the largest objects. Dynamical friction produces a velocity minimum in batches 
that contain the largest fraction of the total mass. 

Fig. 11. — Evolution of the maximum radius with time for models with different initial mass 
(Mq) and initial radius (ro), for low initial eccentricity (e = 10~^). The time scale to reach 
runaway growth decreases with smaller tq and with larger Mq. 

Fig. 12. — Evolution of the maximum radius as in Figure 11, for models with large initial 
eccentricity (e = 10~^). Models with high e require 2-4 times more mass to reach runaway 
growth on time scales similar to that of low e models. 

Fig. 13. — Summary of velocity evolution models for (a) e = 10~^ and (b) e = 10~^. Filled 
circles indicate successful simulations that produce a few Pluto-sized objects and ~ 10^ 
KBOs in 100 Myr or less; open circles indicate simulations that produce no Plutos and too 
few KBOs in 100 Myr or less); filled circles within a larger open circle indicate partially 
successful simulations that produce a few Pluto-sized objects but too few KBOs in 100 Myr 
or less). 

Fig. 14. — The evolution of Mmaj./(m), as a function of time, (a) constant velocity models, 
and (b) with velocity evolution. The rapid increase in M^axK^^) at the latter stages of 
many simulations indicates runaway growth. 
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Fig. 15. — Evolution of the mass distribution for a constant cross-section, Aij — ac- The sohd 
hnes plot analytic results for four values of rj; the symbols indicate results of the numerical 
simulations. 

Fig. 16. — Evolution of the mass distribution for Aij = f3c{Mi + Mj). The solid lines plot 
analytic results for four values of rj; the symbols indicate results of the numerical simulations. 

Fig. 17. — Evolution of the mass distribution for Aij = ■jcMiMj. The solid lines plot analytic 
results for four values of rj; the symbols indicate results of the numerical simulations. 

Fig. 18. — Runaway growth for Aij — ^ycMiMj. (a) Evolution of the residual mass 
distribution for four values of 77 > 1. The simulations lag the analytic model for 77 f=i 1 
and then follow it closely for larger 1]. (b) Evolution of the mass of the runaway body for 
the simulation (symbols) and the analytic model (solid curve) as a function of rj. 
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